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Abstract 

Mixed atomistic and continuum methods offer the possibility of carrying out simulations 
of material properties at both larger length scales and longer times than direct atomistic 
calculations. The quasi-continuum method links atomistic and continuum models through 
the device of the finite element method which permits a reduction of the full set of atom- 
istic degrees of freedom. The present paper gives a full description of the quasicontinuum 
method, with special reference to the ways in which the method may be used to model 
crystals with more than a single grain. The formulation is validated in terms of a series of 
calculations on grain boundary structure and energetics. The method is then illustrated in 
terms of the motion of a stepped twin boundary where a critical stress for the boundary mo- 
tion is calculated and nanoindentation into a solid containing a subsurface grain boundary 
to study the interaction of dislocations with grain boundaries. 

1 Introduction 

A longstanding ambition in the modeling of materials has been that of rationalizing and pre- 
dicting the observed mechanical properties of materials on the basis of an understanding of their 
constituent defects. The advent of ever more powerful computers alternatives has ushered in 
the possibility of carrying out rudimentary calculations of this type directly on the basis of full 
atomistic simulations. However, an alternative class of models has sought to exploit atomistic 
insights without abandoning altogether the powerful resources that are associated with contin- 
uum theories. One such approach is the quasicontinuum method (Tadmor, Ortiz and Phillips, 
1996) which links the kinematic constraints, and attendant degree of freedom reduction offered 
by the finite element method, to total energies predicated entirely on atomistic analysis. 
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The push to develop models of defect interactions has come from experimental observations 
on ever-smaller length scales. Recent micromechanical observations routinely explore problems 
in which small numbers of defects are responsible for the mechanical properties of interest (see, 
for example, Gerberich, Nelson, Lilleodden, Anderson and Wyrobek, 1996). The experimental 
advent of nanomechanics as ushered in by a host of high resolution microscopies such as high- 
resolution TEM and atomic scale resolution surface probes such as the STM and the AFM 
has led to theoretical demands as well. One of the theoretical responses to this challenge has 
been the attempt to build simulation tools that allow for the analysis of multiple length scales 
simultaneously. Before turning to the quasicontinuum method itself, we mention two examples 
that are antecedents to the approach advocated here. 

The broad class of models known as cohesive zone models have as their aim the incorporation 
of constitutive nonlinearity to account either for the "core" material at a crack tip (Barenblatt 
1962) or at a dislocation core (Peierls 1940). One of the significant outcomes of these calculations 
that is especially noteworthy for the present purpose is the fact that the incorporation of the 
constitutive nonlinearity introduced by the cohesive zone eliminates (or at least ameliorates) 
the singularities that are inherited from an analysis predicated purely on the basis of linear 
elasticity. A recent advance in cohesive zone technology has been the exploitation of cohesive 
zone elastic potentials calculated using density functional calculations (Xu, Argon and Ortiz, 
1995), in the study of dislocation nucleation near a crack tip. As will be more evident below, the 
quasicontinuum method takes its cue from the cohesive zone approach in that it is noted that 
constitutive nonlinearity, especially as dictated by an underlying atomistic analysis, that gives 
this approach its power. 

An alternative class of models that is built around the same insights as those associated with 
cohesive zone approaches are those that effect a linkage between two different spatial regions, one 
of which explicitly treats each and every atomic degree of freedom with its requisite nonlinearity 
and a second of which is treated either by traditional linear elasticity or its discrete analog 
(Kohlhoff, Gumbsch and Fischmeister (1991), Thomson, Zhou, Carlsson and Tewary (1992)). 
Such methods necessitate the management of boundary conditions that permit the matching of 
the two regions. However, the logic that stands behind such methods is nearly identical to that 
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advocated here. The assertion is that in the immediate vicinity of the defect, it is essential to have 
a sufficient level of resolution as to capture the sometimes complex local atomic rearrangements 
that characterize them. On the other hand, the contention is that far away from the defect, 
there is no reason to doubt the efficacy of continuum mechanics. In their present incarnations, 
these mixed methods have the difficulty in that each problem must be treated on a case by case 
basis with the boundary between the fully atomistic region and the rest of space designed a 
priori. The quasicontinuum method to be described below is founded on the basis of adaptive 
strategies in which, in the course of the calculation, the region of full atomistic resolution can 
vary in response to the evolution of both loading and defects. 

The present paper is written with two clear objectives. First, it is our aim to generalize the 
earlier statements of the quasicontinuum method so as to allow for the treatment of problems 
involving more than one grain. It is well known that a range of phenomenology in the mechanics 
of materials including Hall-Petch type phenomena, grain boundary sliding and Nabarro-Herring 
creep, trace their existence to the presence of grain boundaries. The formulation as described 
in an earlier paper (Tadmor, Ortiz and Phillips, 1996) was noncommittal with respect to the 
question of how one might incorporate grain boundaries as a synthetic element of the material's 
microstructure and it is a key aim of the present work to remove that ambiguity. 

The second objective of the present paper is to provide the logical foundation for the gener- 
alization of the method to a number of different contexts. Work to further extend the quasicon- 
tinuum method is presently in progress in a number of different areas, all of which rest upon the 
foundations laid here. One of the key extensions is to make the method fully three dimensional, 
which allows for the realistic simulation of phenomena such as dislocation junction formation 
and dislocation nucleation. An additional area of development concerns the need to explicitly 
evaluate the dynamical trajectories undertaken in a given process at finite temperature. In par- 
ticular, we see the quasicontinuum method as a possible meeting point for molecular dynamics 
and continuum thermodynamics. Finally, to allow for the treatment of alloys and even elemen- 
tal materials such as silicon, the method has required generalization to situations in which the 
crystallography is characterized by the presence of more than one atom per unit cell. Although 
these topics will not be described in the present paper explicitly, preliminary efforts have been 
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forged in all of these directions, all of which have been built around the ideas that are set forth 
here. 

The organization of the remainder of the paper is as follows. Section 2 presents an overview 
of the quasicontinuum method. The detailed description of the individual components of the 
methodology may be found in Section 3. This will be followed in section 4 by a discussion of 
validation of the method, with full atomistic calculations serving as the benchmark for success 
with particular reference to the use of quasicontinuum calculations for the study of grain bound- 
ary structure and energetics. Sections 5 and 6 complete the description with applications that 
exploit the capacity of the method to examine deformation problems involving more than a single 
grain. 

2 Overview 

In this section, we undertake a description of the quasicontinuum method, with special attention 
drawn to the subtleties that attended the generalization to problems involving more than one 
grain. As discussed in the introduction, one of the fundamental precepts around which the 
quasicontinuum method was built was the idea that direct atomistic simulation has both strengths 
and weaknesses. The clear advantage of the atomistic perspective is that such calculations are 
capable of providing the requisite resolution to account for the highly nonlinear and sometimes 
counterintuitive atomic arrangements that are found in the defect core. Indeed, it has been 
found that in some circumstances such atomic level details are the source of known macroscopic 
anomalies in the material's behavior (Christian 1983). On the other hand, the weakness of the 
atomistic approach is the huge number of redundant degrees of freedom away from such defects. 
The quasicontinuum method attempts to incorporate both of these insights by allowing for full 
atomic scale resolution near defect cores while exploiting a coarser description further away. 

We begin with an overview of the conceptual elements of the quasicontinuum method as 
described in Tadmor, Ortiz and Phillips (1996). The key idea is that of kinematic slavery in 
which by virtue of the finite element method (FEM), the positions of the majority of atoms are 
entirely constrained and determined only by the displacements of the nodes tied to the element 
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of which they are a part. Once the geometric disposition of the body is established, the problem 
becomes one of determining the total energy. From a traditional continuum mechanics viewpoint, 
such as that offered by the linear theory of elasticity, the total energy may be computed on the 
basis of a phenomenological constitutive law such as Hooke's law. The aim in Tadmor, Ortiz, 
Phillips (1996), by way of contrast, has been to use atomistic calculations to inform the energetic 
statement of the continuum mechanics variational principle. This step offers the constitutive 
nonlinearity alluded to above, and allows naturally for the emergence of lattice defects such as 
dislocations. 

We now give a formal description of the formulation of the quasicontinuum method generalized 
to account for the presence of multiple grains, restricting attention to those problems in which 
the undeformed state of the body is polycrystalline. We take the view that the body whose 
disposition is of interest should be thought of as a collection of a possibly huge number, N, of 
atoms. In fig. |l|, we show the body which is imagined to be built up of a variety of different 
grains with Bravais lattice vectors schematically indicated. The presence of a crystalline reference 
configuration is exploited in the sense that for many regions of the crystal, it is unnecessary to 
save lists of atomic positions since they can be generated as needed by exploiting the crystalline 
reference state. A given atom in the reference configuration is specified by a triplet of integers 
1 = {h,h,h), and the grain to which it belongs. The position of the atom in the reference 
configuration is then given as, 

x(o = EW + RM > (i) 

a=l 

where is the a th Bravais lattice vector associated with grain and R M is the position of a 
reference atom in grain which serves as the origin for the atoms in grain G^. 

Once the atomic positions have been given, from the standpoint of a strictly atomistic per- 
spective, the total energy is given by the function 

Etot = E exact (-xi, x 2 , x 3 , xtv) = E exact ({-Xi}) , (2) 

where Xj is the deformed position of atom i. We adopt the convention that capital letters refer 
to the undeformed configuration while lower case letters refer to the deformed configuration. 
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The energy function in eq. (§) depends explicitly upon each and every microscopic degree of 
freedom, and as it stands becomes intractable once the number of atoms exceeds one's current 
computational capacity. The problem of determining the minimum of the potential energy is in 
the context noted above nothing more than a statement of conventional lattice statics. We will 
now proceed to construct the formulation of the method in a step-by-step fashion. We begin with 
the introduction of kinematic constraints which have the effect of reducing the number of degrees 
of freedom being accounted for by selecting a small subset of R atoms from the total set of N 
atoms. These atoms serve as "representative atoms" , and remain the only unconstrained degrees 
of freedom in the problem. A finite element mesh with nodes corresponding to the positions 
of the representative atoms is then defined. By virtue of finite element interpolation we may 
compute the total energy from the equation given above, but now with a substantial fraction 
of the atoms participating geometrically as nothing more than drones since their positions are 
entirely determined by the displacements of the adjacent nodes. Quantitatively, if a is an index 
over the representative atoms, then the interpolated position x* n * of any other atom i may be 
obtained by 

xf*=I>a(X0x a , (3) 

a 

where iV Q (Xj) is the finite element shape function centered around the representative atom a 
(which is also a FEM node) evaluated at the undeformed position Xj of atom i. In particular, 
we may write the total energy as 

Etot — Eexacti^X 1 ^ , X^™*, ...,X^*) = E exact ( { X*™' } ) . (4) 

At this stage not much has been gained since the computation of the total energy is still 
predicated upon a knowledge of all of the atomic positions, though now many of the atomic 
positions are constrained. We now make the additional assumption that the energy may be 
written in a form that is additively decomposed, such that 

N 

E tot = Y, E n ( 5 ) 
i=i 

which presupposes the existence of well-defined site energies Ei, and is typical in many current 
atomistic formulations such as the embedded-atom method (EAM). The summation runs over 
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all atoms in the solid. Because of this, the sum written above remains intractable in the sense 
that if we interest ourselves in computing the total energy, we are still obliged to visit each 
and every atom. The spirit of the problem that we are faced with is now identical to that of 
numerical quadrature, and what we require at this point is a scheme for approximating the sum 
given above by summing only over the representative atoms with appropriate weights selected 
so as to account for differences in element size and environment. In particular, we desire 

R 

Etot ~ Yl n « E <*- ( 6 ) 

a=l 

The crucial idea embodied in this equation surrounds the selection of some set of representative 
atoms, each of which are intended to characterize the energetics of some spatial neighborhood 
within the body as indicated by the weight n a . As yet, the statement of the problem is incomplete 
in that we have not yet specified how to determine the summation weights, n a . We treat 
the problem of the determination of n a in a manner analogous to determination of quadrature 
weights in the approximate computation of definite integrals. In the present context the goal 
is to approximate a finite sum ("definite integral" on the lattice) by an appropriately chosen 
quadrature rule where the quadrature points are the sites of the representative atoms. Physically 
the quantity n a may be interpreted as the "number of atoms represented" by the representative 
atom a. The quadrature rule (eq. (El)) is designed such that in the limit in which the finite element 
mesh is refined all the way down to the atomic scale (a limit we denote as fully refined) each 
and every atomistic degree of freedom is accounted for, and the quadrature weights are unity 
(each representative atom represents only itself). On the other hand, in the far field regions 
where the fields are slowly varying, the quadrature weights reflect the volume of space (which is 
now proportional to the number of atoms) that is associated with the representative atom. The 



details of this procedure may be found in section |37[ 



The description given above describes the essence of the formulation as it is currently prac- 
ticed. We now describe an additional energetic approximation that simplifies the energy calcu- 
lations and also makes it possible to formulate boundary conditions which mimic those expected 
in an elastic continuum. The essential idea is motivated by fig. 0, which depicts the immediate 
neighborhood of a dislocation core. In particular, for this Lomer dislocation we note the char- 
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acteristic geometric signature of the core, namely, the pentagonal group of atoms in the core 
region. We now focus our attention on the environments of two of the atoms in this figure, one 
(labeled A) in the immediate core region, and the other (labeled B) in the far field of the defect. 
It is evident that the environment of atom A is nonuniform, and that each of the atoms in that 
neighborhood experiences a distinctly different environment. By way of contrast, atom B has an 
environment that may be thought of as emerging from a uniform deformation, and each of the 
atoms in its vicinity sees a nearly identical geometry. 

As a result of these geometric insights, we have found it convenient to compute the energy 
E a from an atomistic perspective in two different ways, depending upon the nature of the atomic 
environment of the representative atom a. Far from the defect core, we exploit the fact that the 
atomic environments are nearly uniform by making a local calculation of the energy in which it 
is assumed that the state of deformation is homogeneous and can be well-characterized by the 
local deformation gradient F. To compute the total energy of such atoms, the Bravais lattice 
vectors of the deformed configuration b a are obtained from those in the reference configuration 
B a via b a = FB a . Once the Bravais lattice vectors are specified, this reduces the computation 
of the energy to a standard exercise in the practice of lattice statics. 

On the other hand, in regions that suffer a state of deformation that is nonuniform, such as the 
core region around atom A in fig. |2], the energy is computed by building a crystallite that reflects 
the deformed neighborhood from the interpolated displacement fields. The atomic positions of 
each and every atom are given exclusively by x = X + u(X), where the displacement field u 
is determined from finite element interpolation. This ensures that a fully nonlocal atomistic 
calculation is performed in regions of rapidly varying F. An automatic criterion for determining 
whether to use the local or nonlocal rule to compute a representative atom's energy based on the 



variation of deformation gradient in its vicinity will be presented in section |3.2| . The distinction 
between local and nonlocal environments has the unfortunate side effect of introducing small 
spurious forces we refer to as "ghost" forces at the interface between the local and nonlocal 



regions. A correction for this problem is also discussed in section |3T2 . 

Once the total energy has been computed and we have settled both the kinematic and ener- 
getic bookkeeping, we are in a position to determine the energy minimizing displacement fields. 
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As will be discussed below, there are a number of technical issues that surround the use of either 
conjugate gradient or Newton-Raphson techniques. Both of these techniques are predicated upon 
a knowledge of various derivatives of the total energy with respect to nodal displacements and 



are explained in section 3.3 



As noted in the introduction, one of the design criteria in the formulation of the method 
was that of having an adaptive capability that allowed for the targeting of particular regions 
for refinement in response to the emergence of rapidly varying displacement fields. For example, 
when simulating nanoindentation, the indentation process leads to the nucleation and subsequent 
propagation of dislocations into the bulk of the crystal. To capture the presence of the slip that 
is tied to these dislocations, it is necessary that the slip plane be refined all the way to the atomic 
scale. The adaption scheme to be described in section p.4| allows for the natural emergence of 
such mesh refinement as an outcome of the deformation history. 

The goal of the present section has been to elucidate the key conceptual elements involved in 
using the quasicontinuum method. However, as was noted above, certain features in the formu- 
lation involve subtleties demanding further attention. In the following sections, we undertake a 
more detailed analysis of some of those issues. 

3 Details of the Methodology 
3.1 Reduced Atomic Representation 

In this section we discuss the selection of the representative atoms, and the construction of an 
expression of the total energy that depends only on the degrees of freedom of the representative 
atoms. We consider the undeformed body to be a crystalline solid, i.e., a collection of N atoms, 
which occupy a region Bq and may be arranged in many grains (see fig. |l|). The undeformed 
position Xj of any atom i is obtained from the coordinate of a reference atom and an associated 
set of Bravais lattice vectors that are specified for each grain as discussed earlier in (eq. ([I])). The 
deformed configuration of the body is described by a displacement function u which depends on 
X and the deformed position of any atom i can be obtained as 

Xj = Xj + Uj, (7) 
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where = u(X.j). 

On loading the solid, the equilibrium configuration of the body is defined by the set of 
displacements Uj which minimizes the potential energy function 

N 

Il(u) = £fo t (ui, ...,Ujv) - X] f i • U i> ( 8 ) 

i=l 

where E tot is the total energy of the system obtained from an atomistic formulation, fj is the 
external force acting on the atom i and u satisfies the essential boundary conditions of the 
problem. This is the well-known method of lattice statics (LS). Now, as stated earlier, we 
assume that E tot can be decomposed as a sum over the energies of individual atoms E^, i.e., 

N 

E tot (ut, uat) = Ei(ut, ujv), (9) 

i=i 

and eq. (||) becomes 

N N 

n(u) = E i("i, •••> "n) - £ f i ■ u - ( 10 ) 

i=i i=i 

Many of the conventional atomistic formulations (such as the embedded-atom method) admit 
such a decomposition, although it is not admissible in the more sophisticated density functional 
approach. For example, the embedded-atom method (Daw, Baskes 1986) provides that for a 
homonuclear material the energy at site % is given by 

E l = WHn 1 ) + f(p l ), (ll) 

Z 3 

where is the distance from atom % to the neighbor j, is a pairwise interaction, pi is the 
electron density at the site i and f(p) is the embedding energy. The potentials are assumed 
to have a cutoff radius of r cut , i.e., any atom interacts directly only with those atoms within a 
distance r cut from it. 

The variational principle associated with eq. (|T0|) provides the solution Uj and may lead to 
nonlinear minimization problems with intractable numbers of degrees of freedom. This motivates 
the need to formulate approximation strategies that preserve the essential details of the problem 
while requiring fewer degrees of freedom. The first step in our approximation method is the 
selection of a subset of R atoms (R N) called representative atoms, to describe the kinematics 
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and energetics of the body. To motivate the reasoning underlying this approach, we revisit 
fig. discussed earlier where the atomic structure in the vicinity of a relaxed Lomer dislocation 
in fee aluminum is presented. Two atoms have been highlighted along with their respective 
environments. Atom A lies at the dislocation core itself, while atom B is further away. In 
the region around atom A, the deformation fields are changing rapidly on the scale of atomic 
distances. The non-uniform nature of the deformation near the dislocation core implies that all 
atoms in that region experience completely different environments. At the same time, many of 
the atoms in the vicinity of atom B experience environments nearly identical to that of atom 
B, and thus are nearly energetically equivalent to atom B. This conclusion naturally leads to 
the concept of representative atoms mentioned above. In this case, atom B could well represent 
several of its neighbors, due to their similar environments, while all atoms near A have to be 
chosen as representative atoms. Where the deformation gradients are large and quickly varying 
on the atomic scale more representative atoms will be selected, while further away fewer will be 
explicitly considered. Such a representation is presented in fig. |3]a for the same dislocation core 
structure. We see that near the core region all atoms are represented, while further away, where 
the deformation is more homogeneous, the density of representative atoms is reduced. 

The displacements {u a } of the representative atoms are the relevant degrees of freedom of 
the system. The next task is to construct an approximate energy function Uh that depends only 
on {u Q }. To achieve this we first need a kinematic description of the body, i.e., we need a tool 
to describe the deformed positions of every atom in the body when the displacements of the 
representative atoms are known. This is achieved by the construction of a finite element mesh 
(f2 e ,e = 1...M, where M is the number of elements) with the representative atoms as nodes 
(cf. fig. |3|b). The deformed position of any atom in the model is then obtained by interpolation 
using the finite element shape functions and the displacements of the representative atoms; the 
deformed configuration of the body may thus be completely described. We have chosen to use 
linear triangular finite elements (i.e., the deformation gradient is constant in each element), which 
are generated using the constrained Delaunay triangulation (Sloan 1993). This triangulation 
allows for the convenient treatment of non-convex, multiply connected regions. 

The energy of the body is described by eq. (|6]), and requires a knowledge of n a . As explained 
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above, these quantities may be thought of as quadrature weights for the summation on the 
discrete lattice. We now discuss the computation of n a . Let g be a real valued function on the 
lattice with g(Xj) being its value at the site i. We define 

N 

S(g)=Y,90Q, (12) 

i=i 

as the sum of g. If this sum is to be evaluated using the value of g only at the representative 
atoms, g(X a ), and a quadrature rule, we have 

R 

S h (g) = Y,n a g(X a ). (13) 

a=l 

The values of n a may now be determined by insisting that 

S(gp) = S h {gp) (14) 

for some functions gp {(3 = 1...R) chosen a priori, i.e., by insisting that the quadrature rule 
should sum the functions gp exactly. There are many possible choices of gp, some of which we 
discuss below. 

1. Voronoi characteristic functions: The function gp is chosen to be the characteristic function 
of the Voronoi cell Vp (Okabe, Boots and Sugihara, 1992), associated with representative 
atom p, 

X V p(*i) = 1, VXiG^g 

= 0, otherwise. (15) 

Setting gp = Xp and applying the relation in eq. ([14]) we obtain rip to be 

np = S h (x V p)=S( X V p)) = j:x V p(^). (16) 

i 

It is easily seen that n a admits a simple interpretation as the total number of atoms in the 
Voronoi cell of representative atom a. 

2. Patch characteristic functions: We define a patch (Pp) surrounding a representative atom/node 
(3 as the polygon constructed by joining the centroids and midpoints of the edges of the 
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elements incident on the representative atom. The function gp is chosen to be the charac- 
teristic function of the patch Pp and it is now immediate that 

N 



n p = S(x P ,)=Y,X P ,(^) (17) 



i=l 



and again np may be interpreted as the number of atoms contained in the patch Pp. 

3. FE shape functions: Here we choose the function gp to be the finite element shape function 
Np, and in a manner similar to the above two cases it follows that 

n p = S(N p )= y £N p {X i ). (18) 



Some remarks are in order. First, note that in each case, the functions gp are chosen so that 
their value is unity at the representative atom f3 and vanishes at all other representative atoms. 
Second, the first two alternatives require the explicit construction of a tessellation (either Voronoi 
cells or the patches). Third, all the alternatives provide the quadrature weight to be unity at a 
representative atom situated in a fully-refined region. It has been found that the results of the 
calculations are insensitive to the choice of the different schemes. 

Armed with the values of n a and eq. @ the approximate energy function may be written 

as 

R R 
U h (u) = Wa-So(Ul, U R ) - n oS-cx ■ U Q , (19) 
a—1 a=l 

where f is the average force on representative atom a, and the subscript h refers to the approxi- 
mation introduced by the finite element partitioning. In the event that one chooses all the atoms 
in the model the energy in eq. (|19D reduces to the exact function of eq. fllCf ) and one recovers 
lattice statics. 

3.2 Computation of Representative Atom Energies 

The energy of any representative atom may be computed by creating a list of its neighbors (we call 
such a list the representative crystallite) , and obtaining the deformed positions of these neighbor 
atoms using the finite element interpolation. This approach will require the explicit neighbor 
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list computation of each of the representative atoms and proves to be very time consuming. A 
more efficient strategy can be formulated as follows. Consider the atom B in fig. |2j. This atom 
experiences only a slowly varying deformation, and its energy can be well approximated by that 
computed using the local deformation gradients and the Cauchy-Born rule (Ericksen 1984 and 
references therein) which states that the atoms in a deformed crystal will move to positions 
dictated by the existing gradients of displacements. 

On the other hand, in the region around atom A, the deformation fields are changing rapidly 
on the scale of atomic distances. This observation suggests a division of the representative atoms 
into two classes (a) Nonlocal atoms whose energies are computed by an explicit consideration of 
all its neighbors and (b) Local atoms whose energies are computed from the local deformation 
gradients using the Cauchy-Born rule. The former type of representative atoms are essential 
to capture the atomistic nature of defect cores and interfaces while the latter represents the 
continuum limit and is used in regions of the solid undergoing only a near homogeneous de- 
formation. It is emphasized that the local formulation is an approximation that provides for 
efficient computation. 

We now discuss the criterion that decides if a given representative atom is to be treated 
as local or nonlocal. The representative crystallite of a representative atom experiences various 
deformation gradients arising from the different elements that surround it. A state of deformation 
near a representative atom is near homogeneous if the deformation gradients that it senses from 
the different elements are nearly equal. This can be characterized by the inequality 

max|A£-A*| < e L , (20) 

a,b;k 

where a and b run over all elements that are within some radius r s of the representative atom 
discussed presently, \ e k is the k eigenvalue of the right stretch tensor U e (see for example, 
Chadwick 1976) obtained from the deformation gradient F e in element e, and is an empirically 
selected constant. This criterion implies that the error in the computation of the energy using 
the local deformation gradients when compared with a fully atomistic calculation of the energy 
is within a specified tolerance. Thus all atoms for which eq. (|20| ) is satisfied are treated as local 
atoms while for the remaining atoms the nonlocal rule is used to compute the energy. In addition, 
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any atom which is within r s of a surface or interface of interest is made nonlocal. 

Using the above criterion, the total number of nonlocal atoms in the model is determined 
by two prescribed parameters, the tolerance on the eigenvalues, e^, and the range of nonlocal 
influence, r s . For correct surface and interfacial energy, atoms within r cut of the interface must be 
nonlocal (i.e. r s = r cut ). For correct forces and stiffness (promising a correct relaxed configura- 
tion) each of these atoms must be embedded in a further r cut radius of nonlocal atoms, bringing 
the total necessary nonlocal padding to r s = 2r cut . A lesser padding will result in less demand- 
ing computational models with fewer degrees of freedom at the expense of more approximate 
interfacial and defect core structures. 

The total energy (eq. (|l9|) ) can now be separated into its local and nonlocal contributions, 

Rl Rnl R 

n h (u) = £ n a E l ° c {¥ Xl F M ) + npEpim, u R ) - £ nj a ■ u a , (21) 

a=l j3=l a=l 

where there are Rl local atoms and Rnl nonlocal atoms (such that Rl + Rnl = R)- 

The procedure used in the calculation of energies of the local atoms from the local gradients 
of deformation is taken up presently. Consider representative atom a experiencing near homo- 
geneous deformation. If n e a denotes S(x{Q e ) • ipa), which may be interpreted as the number of 
atoms associated with the element fl e represented by the representative atom a, then the energy 
E l ° c is approximated by 

M M 

n a E l a oc = J2<£(F e ), n a = Y,<, (22) 

e=l e=l 

where £(F) is the energy of a single atom experiencing a homogeneous deformation given by the 
deformation gradient tensor F. 

The nonlocal energy Ep(ui, u^) is computed as it would be in a standard atomistic cal- 
culation. The energy of atom (3 is a function of the positions of all atoms falling inside its 
cutoff sphere for the given deformation. We assume there are mp such neighbors which occupy 
positions X^,j = 1,2, ...,mp in the undeformed configuration. The deformed positions of these 
atoms relative to atom f3 will be, 

rj = X£ + 14 - X, - (23) 
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where and are the coordinates and displacements at atom j3 and is the displacement 
interpolated to the site of neighbor j of atom (3 and is given by, 

uj=f;u a JV a (Xj), (24) 

a=l 

where u Q is the displacement at representative atom/node a, N a is the associated finite element 
interpolation function. Substituting eq. ( p2|) into eq. ( pTj ) and specifically accounting for the 
dependence of Ep on we have, 

«l M Rnl R 

n ft (u) = E E <£( F e ) + E n/j^(rj, i7) - E n«A • u, (25) 

a=le=l 3=1 a=l 

The sums in the first expression of eq. ( p5|) can be reversed so that, 

M Rnl R 

U h (u) = E fe£(F e ) + E ^( r J> -> r 7) - E • U «> (26) 

e=l /3=1 a=l 

where 

Rl 

= E < (27) 

is the total number of atoms (associated with local representative atoms) falling in element e. 

Although eq. ( p6|) is a complete description of the approximate energy function, and represents 
a mathematically consistent formulation, it leads to noisy solutions in the presence of nonlocal 
representative atoms with weights exceeding unity. Expressed differently, solutions are smoother 
when nonlocal representative atoms are present only in the fully refined regions. The cause of 
this noise or error may be traced to the non-uniformity in the constrained kinematics due to a 
non uniform mesh (this effect has been noted by Tadmor, Ortiz and Phillips (1996) and detailed 
in Tadmor (1996)). The remedy for this problem is to fully refine the mesh around atoms that 
are treated nonlocally, and as a result the weights associated with the nonlocal atoms will be 
unity. The resulting approximate energy function eq. (EST) reduces to 



M Rnl R 

n h (u) = E ^(F e ) + E M4> r 7) - E nJa ■ u«, (28) 

e=l 13=1 a=l 

In the type of configurations we have studied, nonlocal atoms tend to appear in groups we 
refer to as clumps, surrounded by local atoms. The local and nonlocal formulations are not 
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completely compatible and, although there is no seam in the formulation of the energy function, 
non-physical forces arise on atoms in transition zones between local and nonlocal regions, even 
when the crystal is undeformed. To illustrate this point, consider a fully refined mesh as shown in 
fig. |j. The atoms represented by open circles are local atoms whereas the dark ones are nonlocal. 
The circle drawn around nonlocal atom A represents its cutoff sphere. On the one hand, as local 
atom B lies inside the sphere, its motion affects A and dE^ L /<9ub 7^ 0, which corresponds to 
a force acting on B due to A. On the other hand, the energy of B is computed locally and 
depends only on atoms which share a common element with B (the shaded elements in fig. f|). 
Therefore, dE^/du^ = and the motion of A does not lead to forces on B. These imbalances 
result in non-physical forces on A and B, which we call "ghost" forces. Their order of magnitude 
is 0.1 eV/A and can lead to energy relaxation of order 0.005 eV. The principal reason for this 
erroneous unbalanced forces is the fact that our procedure has focussed on approximating the 
energy and not the forces. 

In order to avoid the ghost forces, we have to correct the forces acting on atoms in the 
transition zones. This is achieved by demanding that forces acting on any atom be computed 
using only the formulation which corresponds to its status, as if the atom was in a fully local or 
nonlocal region. For example, the nonlocal term dE^ L /<9ub is not added to the forces acting on 
local atom B and a term dE^ L /8ua is computed for atom A. Ghost forces don't derive from 
a potential and as a result they are not symmetrical i. e., the ghost force on B due to A is not 
equal to that due to A on B. Thus they cannot be corrected in the energy directly. They are 
instead computed each time the status of the representative atoms are updated and are then 
held constant until next update. These dead loads f G can be incorporated as corrections to the 
energy function: 

n' /l (u) = n,(u)-5:f Q G -u Q . (29) 

a 

The ghost forces are a function of atomic positions and thus new ghost forces arise, once 
atoms relax. Their norm is linked to the size of relaxation in transition regions. For example, 
in the case of the relaxation of a twin boundary in aluminum (see section f|), the difference in 
ghost forces before and after relaxation is 3 x 10~ 6 eV/A corresponding to an energy relaxation 
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of order 10~ 7 eV. More generally, we observe that the magnitude of the ghost forces is decreased 
only by a factor of the order of ten when the dead load approximation described here is applied. 
Nevertheless, this procedure improves the accuracy of the solutions in transition zones. 
The complete energy function including the approximate ghost force correction is then, 

M Rnl R 

II'Ju) = 5> e £(F e ) + £ Efi {$, ...,r7) - £ n a {f* + f Q G ) • u a . (30) 

e=l /3=1 o=l 

3.3 Energy Minimization 

We are interested in obtaining equilibrium configurations of the solid. We thus invoke the 
principle of minimum potential energy which states that a system will be at equilibrium when 
its potential energy is minimum. The potential energy of our reduced atomic system is given 
above in eq. (|3~0f). To minimize this energy we have used both conjugate gradient algorithms 
(Papadrakakis and Ghionis 1986) and quasi-Newton algorithms (Denis and Schnabel 1983). The 
conjugate gradient algorithm requires computation of the gradient of eq. (|30f) with respect to 
the representative atom displacements (i.e., the system degrees of freedom). The quasi-Newton 



method requires in addition also the Hessian of eq. (|30|) , i.e. the second gradient with respect to 
displacements. These will be evaluated in this section for generic interatomic interactions that 
can be written in the form given by eq. (|5|). 

The gradient of the potential energy (also referred to as the out-of-balance force vector) is 
given by, 



Ej or f3 



njf a + ff ), (31) 

/ f ii n 

3=1 

where P = dS/dF is the first Piola-Kirchhoff stress tensor and ip^ = —dEp/dr^ is the force on 
atom j3 due to its neighbor j. The deformation gradient F e can be expressed in terms of the 
nodal displacements and finite element interpolation functions as, 

R 

F e = I + ]Tu a VoiV a (X e ), (32) 

a=l 

where Vo = <9/9X is the material deformation gradient, I is the identity matrix, and for the 
constant strain elements we use, X e is the element centroid. We then have, 

dF 

^=V iV a (X e ), (33) 



and similarly by using eq. (|23|). 
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(jV (Xj) - 8J) I. 



(34) 



Substituting eqs. (|33|)-(|3~4]) into eq. fl3T|) and rearranging we have, 



^TT' M R NL 

h - ]>> e P(F e )V iV a (X e ) - ]T 



du 



e=l 



(3=1 



rap 

E<^( x £) 



(35) 



3=1 



where it is understood that the third sum is zero when atom a is local. In the limit where all 
atoms are local the ghost force contribution drops out and eq. fl3"5] ) reduces to, 

9W h \ 



= ]>> e P(F e )VoiV Q (X e ) - nj a . 

loc e=l 



(36) 



We thus have a continuum representation of the boundary value problem with the exception that 
the constitutive law, P(F) = dS/dF, is nonlinear and obtained from an atomistic calculation. 

In the fully-refined limit where all atoms are represented and nonlocal, the ghost forces 
similarly drop out, and so does the term with the local contributions. We can separate the first 
nonlocal sum into two parts, 



nonloc 



Rnl 

E 

0=1 



E^a(X 

3=1 



3=1 3=1 



(37) 



Here the shape function acts as a Kronecker delta function, thus iV Q (X^) is equal to one when 
neighbor j of atom (3 is atom a. We thus see that the second sum drops out since N a (K J a ) is 
always zero (atom a cannot be a neighbor of itself), so eq. ([J7]) reduces to 



dm 



nonloc 



Rnl 

-E 

0=1 



nip 

3=1 



+ E^"fa 

3 = 1 



(38 



The first term is the sum of forces exerted on all other atoms by atom a. The second term is 
the sum of forces exerted on atom a by all of its neighbors. The third term is the external force 
acting on atom a. 

For conjugate gradient minimization the expression in eq. (^) is all that is needed. The 
algorithm constructs a series of conjugate search directions from the current gradient and previous 
search directions and proceeds to minimize along each direction using a line search routine (see 



19 



Papadrakakis and Ghionis (1986) for more details). An alternative approach is to iteratively 
solve dH' h /du a = by substituting, 



dm 



dm 



i+i 



+ E 



d 2 W 



du a dur. 



5u 



i+i 



0. 



(39) 



where / is an iteration counter. The Hessian or second gradient expression is the stiffness matrix 
K a/ g and is obtained by differentiating eq. (|35|). Following a similar procedure to that used for 
obtaining eq. (^) itself we have, 



M 



Rnl 



K 



a/3 



$> e C(F e )V iV a (X e )VoiV>(X e ) + J2 



e=l 



EE^P^K 



7=1 Lfc=l 1=1 

rrifj rap 



E E * k MK) -EE 4 lN ^) + ^ E E < 



(40) 



fc=l 1=1 k=l 1=1 k=l 1=1 

where C = d 2 S/dF 2 is the Lagrangian tangent stiffness tensor and k^ 1 = <9 2 i^/<9rg<9r^ is an 
atomic level stiffness. 

The solution then proceeds iteratively by solving eq. (|39|) as stated or in reduced notation, 



£Ky U J +1 + Fi = 0, 



(41) 

where F^ = (dU' h /du a )\i is the out-of-balance force vector defined above in eq. (|35|). In the 



Newton-Raphson method the displacements at each iteration are updated by, 



u 



ul + 5u 



i+i 



(42) 



while for a quasi-Newton solver a line search minimization is done along the search directions 
given by 5u^, +1 . The procedure continues until ||F^|| is reduced sufficiently for all a. 

3.4 Automatic Adaption 

The realization that much of the computation in straightforward atomistic simulation is wasted 
due to the sufficiency of local continuum approximations far from defects is not new. A number 
of mixed continuum and atomistic models have been proposed in recent years to capitalize on this 
feature (some were referenced in the introduction and others can be found in Tadmor 1996). The 
standard approach in these models is to a priori identify the atomistic and continuum regions and 
tie them together with some appropriate boundary conditions. In addition to the disadvantage of 



20 



introducing artificial numerical interfaces into the problem, a further drawback of these models 
is their inability to adapt to changes in loading and an evolving state of deformation. Take for 
example the problem of nanoindentation. As the loading progresses and dislocations are emitted 
under the indenter, the computational model must be able to adapt and change in accordance 
with these new circumstances. 

In the current formulation we tie the need for automatic adaption to an estimate of the error 
introduced by the reduction of degrees of freedom. It is then possible to identify regions where 
this error estimator is high, and subsequently add degrees of freedom in these regions. The result 
is an automatic adaption scheme analogous to adaptive remeshing in finite elements. 

To include such an adaption procedure in the method, we appeal to finite element literature, 
where error estimators and automatic mesh refinement have been subjects of extensive research. 
Recall that our collection of representative atoms are also nodes on a finite element mesh of con- 
stant strain triangles. Thus, we use the error estimator first introduced by Zienkiewicz and Zhu 
(1987) in terms of stresses and later modified by Belytschko and Tabbara (1993) to estimate er- 
rors in the strain fields. In our case, the deformation gradient, F is already needed for computing 
atomic energies, forces and stiffness, and therefore it is convenient to write the Zienkiewicz-Zhu 
error estimator directly in terms of the deformation gradient. Thus, we define the discretization 
error in element e as 

-, 1/2 

/ (F — F e ) T (F — F e )d£l , (43) 

where F e is the finite element solution for the deformation gradient in element e, and F is the 
/^-projection of the finite element solution for F, given by 

F = Nf . (44) 

Here, N is the shape function array, and f is the array of nodal values of the projected deformation 
gradient F. Because the deformation gradient is constant within each constant strain element, 
the nodal values f are simply computed by averaging the deformation gradients over all of the 



elements in contact with the node of interest. The integral in equation eq. (|43l) can be computed 
quickly and accurately using a three-point Gaussian quadrature rule. Elements for which the 
error e e is greater than some prescribed error tolerance are targeted for refinement. 
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Refinement then proceeds by adding three new representative atoms at the atomic sites closest 
to the midsides of the targeted elements. Notice that since representative atoms must fall on 
actual atomic sites in the reference volume B , there is a natural lower limit to element size. If 
the nearest atomic sites to the midsides of the elements are the atoms at the element corners, the 
region is fully refined and no new representative atoms are added. Actual examples of evolving 
mesh refinement is given in fig. [5] for the indentation problem. 

In addition to mesh refinement, mesh coarsening is also an important requirement. For exam- 
ple, consider the passage of a dislocation. As the dislocation moves it leaves a trail of fully refined 
mesh in its wake corresponding to previous core positions. Far behind the dislocation the solid is 
undistorted and the high mesh resolution is unnecessary and could be coarsened. To coarsen the 
following algorithm is applied: (1) For each local node/atom the elements surrounding the node 
and the polygon defined by their outer sides is identified. (2) If none of these elements satisfy 
the adaption criterion, remove the current local node and create a new Delaunay triangulation 
of the outer polygon. (3) If none of the new elements satisfy the adaption criterion then the local 
node and all the old elements connected to it are deleted and the new elements are accepted. 
Essentially, the idea is to examine the necessity of each node. To prevent excessive coarsening 
of the mesh far from defects, the nodes corresponding to the initial mesh are usually protected 
from deletion. 

3.5 Putting it all together 

To demonstrate the steps involved in an adaptive quasicontinuum analysis consider the problem 
of nanoindentation depicted in fig. []. Here a rigid rectangular indenter, infinite in the out- 
of-plane direction, is pressed into the free surface of a thin film aluminum single crystal. The 
dimensions and crystallographic orientation are given in the figure. We are interested in applying 
a quasicontinuum analysis to this problem (which is too large to be comfortably tackled by direct 
atomistics) to obtain such data as load versus indentation curves, the criterion for dislocation 
nucleation under the indenter and stress distributions. 

Different boundary conditions are possible to characterize this problem. It is possible to 
model the indenter as well as the film and consider the interactions between indenter atoms and 
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film atoms. However in the interest of simplicity we choose to neglect these effects and model 
the indenter as a rigid displacement boundary condition, i.e. all atoms on the surface under the 
indenter are forced to move down with it. In addition the displacements parallel to the indenter 
face and in the out-of-plane direction can be constrained to mimic perfect stick conditions or 
released for a friction free indenter. The rest of the surface is left free and unconstrained. Far 
from the indenter, symmetry boundary conditions are applied to the model right and left edges 
and the substrate is taken to be rigid with zero displacements at the interface. 

We need to select an initial set of representative atoms. A fully-refined mesh in the vicinity of 
the indenter is desired from the start in order to capture surface effects there and have sufficient 
resolution to accommodate the indenter geometry. The details of the initial mesh generation 
may be found in Tadmor (1996). The simulation proceeds as follows: 

[1] Next load step - the indenter is driven another step into the crystal (0.2A in these simu- 
lations) by rigidly displacing the atoms under the indenter downward by the appropriate 
amount. 

[2] Local/nonlocal status computation - the locality criterion defined in eq. ( ^0|) is evaluated 
for each representative atom and its status is determined. Significant preprocessing is done 
at this stage (such as the storage of atom lists and computation of shape functions) to 
speed up the computations (see Tadmor (1996) for details). 

[3] Ghost force evaluation - ghost forces are computed and applied as dead loads. 

[4] Energy minimization - a quasi-Newton solver is used to iteratively minimize the total 
potential energy and to identify the equilibrium configuration of the system subject to the 
new load step boundary conditions. 

[5] Automatic adaption - all elements in the mesh satisfying the adaption criterion are adapted, 
i.e., divided into smaller elements. Since all nodes must occupy atomic sites a natural 
cutoff commensurate with the lattice spacing prevents indefinite adaption. At the same 
time that elements are being checked for adaption they can also be checked for coarsening 
and removed if they are not necessary. 
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[6] If elements have been added (or removed) in the adaption phase the system will no longer 
be in equilibrium (since the system has changed). In this case return to [2] to obtain the 
new relaxed configuration. 

[7] Output - load/displacement data, displacement contours, stress and strain contours, en- 
ergy, atomic structure, etc. 

[8] Proceed to [1] for next load step. 

A series of snapshots of the progressing adaption for this problem were given in fig. [5[ The 
load-displacement curve computed for an embedded atom model of aluminum due to Ercolessi 
and Adams (1992) is given in fig. 0a. The response is initially linear as predicted by elasticity 
theory, until dislocations are nucleated at a critical load. The atomic structure under the in- 
denter after dislocation nucleation is presented in fig. ^b. A far more detailed discussion of this 
simulation and others for different orientations and indenter geometries are presented in Tadmor 
et al. (1997). 

In section |6| nanoindentation in an aluminum bicrystal is discussed in more detail. There the 
nanoindentation is used as a means for generating dislocations, as a "dislocation gun" , in order 
to probe the interaction of dislocations with grain boundaries. 

4 Validation 

As was noted in earlier work, our aim with the quasicontinuum method is to properly recover 
two limiting cases. On the one hand, one aims to restore conventional atomistic simulation in 
the limit that full atomic resolution is adopted everywhere. The benchmark of such a success 
is whether or not the quasicontinuum results for defect cores are in accord with those obtained 
by conventional atomistic simulation. On the other hand, restoration of the continuum limit in 
the event of only long wavelength deformations is revealed in features such as the appropriate 
dispersion relation for long wavelength elastic waves. 

This section contains the validation of the formulation presented in the previous section. We 
compare the results obtained using QC with those obtained from lattice statics (LS) in the cases 
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of a (111) free surface in aluminum, a twin boundary in aluminum, a £5 boundary in gold, 
and a £99 boundary in aluminum. Aluminum was modeled using the embedded-atom potentials 
developed by Ercolessi and Adams (1993) and gold with the Finnis- Sinclair potentials of Ackland, 
Tichy, Vitek, Finnis (1987). In all the cases, representative atoms within a distance of 2r cut from 
the interface are treated using the nonlocal rule. 

The (111) free surface in aluminum was modeled using a block of atoms with dimension 
114A x 65A. Table [1] shows a comparison of the relaxed energies of the atoms in different layers 
with those obtained from direct atomistics. The energies of the first three layers are in good 
agreement with those computed via LS. The relaxation process brings about a separation of 
0.02lA between layers A and B and the layers B and C are closer by 0.003A which are equal to 
those obtained with LS. 

A block of size 118Ax 222A was used to simulate a twin boundary (£ = 3(112)) in aluminum. 
The twin plane was chosen to be the y — plane as shown in fig. §. The comparison of relaxed 
energies obtained from LS and QC are shown in table |2|. The y-displacement (1*2) obtained from 
QC (with and without the ghost force removal algorithm) is compared with that from LS in 
fig. |9|, where the agreement is seen to be excellent when the ghost forces are corrected. The 
errors caused by the ghost forces are also seen in this figure. 

Ackland, Tichy, Vitek and Finnis (1987) have carried out lattice statics simulations of the 
£5(210) boundary in Au. Fig. pH| a shows a comparison of the structure obtained using QC with 
that obtained using LS where the agreement is seen to be excellent. The value of the grain 
boundary energy computed using QC is 670m J/m 2 which is consonant with the 676m J/m 2 
obtained by Ackland, Tichy, Vitek and Finnis (1987). 

A more complex £99(557) boundary was also simulated using QC; the comparison with LS 
solution is shown in fig. [TOlb. The structure obtained here also agrees well with that of Dahmen, 
Hetherington, Okeefe, Westmacott, Mills, Daw and Vitek (1990) who performed a combined 
theoretical/experimental study of this boundary (by visual inspection, a quantitative comparison 
was not made). 

We conclude this section by noting that in all the cases presented above the QC solution of 
the interfacial structure agrees well with those obtained from lattice statics. The implication 
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of this success is that the method has been shown to be a viable alternative to lattice statics 
when simulating grain boundaries and thus may be used for simulations involving interfacial 
deformation. 



5 Interfacial Motion 

The macroscopic plastic behavior of a solid is a cumulative result of the motion of dislocations. 
In addition, grain boundaries also accommodate plastic deformation by processes such as sliding. 
In such cases it is of interest to study inhomogeneities on the grain boundary, i.e., grain boundary 
defects and their interaction with an applied stress. The first example described in this section 
attempts to investigate the effect of stress on a twin boundary with a step. As a second example, 
we describe simulations where we study dislocations interacting with grain boundaries. 

The QC method may be used to study the interaction of interfacial inhomogeneities with an 
external stress. The goals of such simulations are twofold: (a) the determination of the critical 
load required to induce plastic deformation and (b) the elucidation of the mechanism of such 
deformation. The present example is that of a step on a twin boundary (S = 3(111)) in aluminum 
and its interaction with an applied shear stress. The finite element mesh and the associated 
step geometry are shown in fig. |ll[ The step is subjected to a far-field homogeneous shear 
deformation which is effected by the application of kinematic boundary conditions, equivalent 
to the shear stress, on the boundary of the model. Fig. |T2] a shows the relaxed configuration 
of the step in the absence of applied loads. On application of the load, the solid undergoes a 
near homogeneous deformation, and on attainment of a critical stress r*, the stepped boundary 
undergoes an inhomogeneous deformation; the configuration after this event is shown in fig. 0b. 



An examination of fig. [12] reveals the migration of the twin boundary by the nucleation of two 
^■[112] edge dislocations from the corners of the step. The value of the displacement jump 
computed from the simulation corresponds to the Burgers vector of a y[112] dislocation which 
is equal to 1.646 A in the case of aluminum. On nucleation, the dislocations move towards 
the boundary of the simulation cell and are eventually stopped by it. The critical value of 
nondimensional stress r* / fi is found to lie between 0.031 and 0.036 and is found to be insensitive 
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to the cell size chosen for the simulation. Fig. [H| shows a plot of effective load (norm of the 
reaction vector) versus the net global shear strain, where a drop in the effective load is seen at 
the critical strain. This drop in the load can be estimated using a simple linear elastic theory. 
It is seen that the initial part of this curve is a linear function of the global strain, and at these 
strain levels all the strain is elastic. On the nucleation of the dislocations (at a global strain level 
of 0.036), the net elastic strain falls to 0.028. The load corresponding to this strain obtained for 
the linear part of the curve is 3.52 eV/A while the value of the effective load obtained from the 
simulation at a global strain level of 0.036 is 3.72 eV/A. The value obtained from the simulation 
is expected to be higher than that predicted due to the fact that the dislocations are trapped 
near the boundary of the simulation cell, and thus the elastic strain is not reduced to the value 
obtained from the simple analysis. 

It is interesting to contrast the critical stress r* with a typical Peierls stress for a straight 
dislocation. For example, in the case of a screw dislocation in this metal the Peierls stress 
is 0.00068/2 (Shenoy and Phillips 1997), nearly fifty times smaller than the critical stress for 
advancing the twin boundary. As another comparison to set the scale of the stresses determined 
here, the stress to induce motion of the twin boundary can be compared with that to operate a 
Frank- Read source which is a ~ fJib/L, where L is the width of the source (Hull and Bacon 1992). 
In light of this estimate, the stress to induce motion of the twin boundary is of the same order 
as that to operate a Frank- Read source of width « 356 (where b is a typical Burgers vector). 
Although typical Frank- Read sources are larger than 356 and hence operate at even lower stresses, 
the stress found to stimulate motion of the twin boundary is still significantly smaller than the 
ideal shear strength, and is an example of the "lubricating" effect of heterogeneities in the motion 
of extended defects. 

6 Interaction of lattice dislocations with a grain bound- 
ary 

The interaction of dislocations with grain boundaries has been identified as an important fac- 
tor governing the yield and hardening behavior of solids. For example, the dependence of the 
yield stress on the grain size given by the celebrated Hall-Petch relationship (see, for example, 
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Hirth and Lothe (1968)), is explained using a pile-up model which assumes that dislocations are 
stopped by the grain boundary. In this section we illustrate how the QC method can be used to 
build realistic models that address the issue of the interaction of lattice dislocations with grain 
boundaries. For the specific GB we consider, we confirm the hypothesis that a pile-up will indeed 
occur, and that no-slip transmission takes place across the boundary. 

We study the interaction of 4f [110] dislocations with a E = 7(241) symmetric tilt boundary 
in aluminum. Fig. |I4] shows a bicrystal, the top face (between A and B in fig. [14]) of which 



is subject to a kinematic boundary condition that mimics the effects of a rigid indenter. On 
attainment of a critical load, dislocations are nucleated at the point A, and they move towards 
the grain boundary. We investigate the nature of the interaction of these dislocations and the 
grain boundary through the consideration of the following questions: will the dislocation be 
absorbed by the boundary, and if so what is the result of this process? Will the dislocation cause 
a sufficient stress concentration at the boundary so as to result in the nucleation of a dislocation 
in the neighboring grain? 

On application of the load, the bicrystal undergoes some initial elastic deformation and the 
first dislocation is nucleated when the displacement of the top face reaches 14.2 A. This dislocation 



is driven into the boundary and is absorbed without an increase in the load level. Figs. |i~5]a,b 
shows the configuration of the grain boundary immediately before and after this nucleation event. 
It is seen that the dislocation absorption produces a step on the grain boundary. This process 
may be understood based on the DSC lattice by decomposing the Burgers vector into DSC lattice 
vectors (King and Smith, 1980). In our case, we find that 

^[1,1,0] = £[312] + ^[241] (45) 

where |f [312] is the Burgers vector of a grain boundary dislocation parallel to the boundary and 
y-[341] is perpendicular to the boundary plane. A careful examination of fig. |T5| b reveals that 
y|[312] travels along the boundary, and stops on reaching the end of the nonlinear zone. On 
subsequent loading, another pair of Shockley partials are nucleated when the displacement of 
the rigid indenter is 18. 2A, which again does not result in any significant reduction of the load. 
Unlike the first pair, these dislocations are not immediately absorbed by the boundary, and they 
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form a pile-up ahead of the boundary as shown in fig. |15| c. On additional indentation, these 
dislocations are also absorbed by the boundary. The simulation was terminated at this stage. 

The neighboring crystal shows no significant dislocation activity and thus it may be concluded 
that slip is not transmitted into the neighboring grain across the boundary. The absorption of 
the dislocation resulted in a sliding motion of the grain boundary by the passage of a grain 
boundary dislocation and the formation of a step on the grain boundary. The formation of the 
step appears to result in the increased resistance of the boundary to dislocations, as is clear from 
the fact that a significantly higher stress level had to be attained before the absorption of the 
second dislocation. 

It is worth noting the significant computational saving obtained by the use of the QC method 
for this problem. The number of degrees of freedom used in the QC model was about 10 4 
while a complete atomistic model of this problem would have required more than 10 7 de grees of 
freedom. The QC simulation required about 140 hours on a DEC-Alpha work-station while a 
purely atomistic model would have required a parallel supercomputer. 

7 Conclusions 

This paper was set forth with a few main objectives. First, the the quasicontinuum formalism 
as given here was advanced as a basis for considering problems involving multiple grains. The 
viewpoint adopted is that of thinning of degrees of freedom, with regions far away from defect 
cores treated approximately by virtue of finite element interpolation and associated quadrature 
rules for evaluating the discrete sums needed to obtain the total energy. These ideas also serve 
as the basis for the extension of the method to three dimensions and to the incorporation of 
dynamic effects via a finite temperature algorithm. 

The second main objective of the present paper was to validate the method in the context of 
a number of new problems. In particular, we have seen that the method allows for the treatment 
of interfacial structures and the study of deformation processes that involve interfaces. Other 
problems, such as the formation of dislocation junctions and the interactions of cracks and grain 
boundaries which have also been treated using the ideas presented here will be described in 
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forthcoming papers. 
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±Jd\ CI 


LS (eV) 

J-JkJ 1 C V 1 


OC (pV) 


A 




n 3433 


B 


0.0368 


0.0367 


C 


0.0024 


0.0025 


D 


0.0003 


-0.0003 


E 


0.0000 


0.0000 


F 


0.0000 


0.0000 



Table 1: Comparison of energies of atoms as obtained from LS and QC for the (111) free surface. 
Layer A is the outermost layer. 



Layer 


LS (eV) 


QC (eV) 


A 


-0.0051 


-0.0051 


B 


0.0122 


0.0123 


C 


0.0031 


0.0031 


D 


0.0000 


0.0000 


E 


0.0000 


0.0000 


F 


0.0000 


0.0000 



Table 2: Comparison of energies of atoms as obtained from LS and QC for the twin boundary. 
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